Long range Neel order in the triangular Heisenberg model 
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We have studied the Heisenberg model on the triangular lattice using several Quantum Monte 
Carlo (QMC) techniques (up to 144 sites), and exact diagonalization (ED) (up to 36 sites). By 
studying the spin gap as a function of the system size we have obtained a robust evidence for a 
gapless spectrum, confirming the existence of long range Neel order. Our best estimate is that in 
the thermodynamic limit the order parameter m' = 0.41 ± 0.02 is reduced by about 59% from its 
classical value and the ground state energy per site is eo = —0.5458 ± 0.0001 in unit of the exchange 
coupling. We have identified the important ground state correlations at short distance. 
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Historically the antifcrromagnetic spin-1/2 Heisenberg 
model on the triangular lattice was the first proposed 
Hamiltonian for a microscopic realization of a spin liquid 
ground state (GS) @]: 
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where J is the nearest-neighbors antiferromagnetic ex- 
change and the sum runs over spin-1/2 operators. At 
the classical level the minimum energy configuration is 
the well known 120° Neel state. The question whether 
the combined effect of frustration and quantum fluctua- 
tions favors disordered gapped resonating valence bonds 
(RVB) or long range Neel type order is still under debate. 
In fact, there has been a considerable effort to elucidate 
the nature of the GS and the results of numerical [p|^il"| , 
and analytical Jl^pS] works are controversial. From the 
numerical point of view, ED, which is limited to small 
lattice sizes, provides a very important feature |l: the 
spectra of the lowest energy levels order with increasing 
total spin, a reminiscence of the Lieb-Mattis theorem |i~7| ] 
for bipartite lattices, and are consistent with the symme- 
try of the classical order parameter jfj). However, other 
attempts to perform a finite size scaling study of the or- 
der parameter indicate a scenario close to a critical one 
or no magnetic order at all |3|J^]. 

The variational Quantum Monte Carlo (VMC) allows 
to extend the numerical calculations to fairly large sys- 
tem sizes, at the price to make some approximations, 
which are determined by the quality of the variational 
wavefunction (WF). Many WF have been proposed in 
the literature [pPJlO|| and the lowest GS energy estima- 
tion was obtained with the long range ordered type. In 
particular, starting from the classical Neel state, Huse 
and Elser Q introduced important two and three spin 
correlation factors in the WF: 
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where |a;) is an Ising spin configuration specified by as- 
signing the value of Sf for each site and 



tt(z)=T(z)exp [*y(E 5 * 2 -E S 
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represents the three sublattices (say A, B and C) classical 
Neel state in the xy-plane multiplied by the three spin 
term 



T(x) = exp(z/3 ]T li jk SIS]S* k 
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defined by the coefficients jijk — 0, ±1, appropriately 
chosen to preserve the symmetries of the classical Neel 
state, and by an overall factor (3 as discussed in Ref. [[| . 
Since the Hamiltonian is real and commutes with the z- 
component of the total spin, S* ot , a better variational 
WF on a finite size is obtained by taking the real part of 
Eq. (||) projected onto the S* ot = subspace. 

For the two body Jastrow potential v(r) it is also possi- 
ble to work out an explicit Fourier transform v q , based on 
the consistency with linear spin wave (SW) results and a 
careful treatment of the singular modes coming from the 
SU(2) symmetry breaking assumption |l^|nj. This anal- 
ysis gives v q = 1 — v/l + 27 g /l — j q for q ^ and other- 
wise, where 7 9 = [cos (q x ) + 2 cos (q x /2)cos (VSqy/2)] /3 
and the q-momenta are the ones allowed in a finite size 
with A-sites. For a better control of the finite size ef- 
fects we have chosen to work with clusters having all the 
spatial symmetries of the infinite system Q . 

In the square antiferromagnet (AF) the classical part 
by itself determines exactly the phases (signs) of the GS 
in the chosen basis, the so called Marshall sign. For the 
triangular case the exact phases are unknown and the 
classical part is not enough to fix them correctly. There- 
fore, one has to introduce the three-body correlations of 
Eq. Although these do not provide the exact answer, 
they allow to adjust the signs of the WF in a non trivial 
way without changing the underlying classical Neel or- 
der. To this respect it is useful to define an average sign 
of the variational WF relative to the normalized exact 
GS |Vo) as 
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with ip(x) = (x\ip). 

We have compared the variational calculation with the 
exact GS obtained by ED on the N — 36 cluster. For 
completeness we have considered the more general Hamil- 
tonian with exchange easy-plane anisotropy a, ranging 
from the XY case (a = 0) to the standard spin isotropic 
case (a = 1). As shown in Tab. |, in the variational 
approach the most important parameter, particularly for 
a — ► 1, is the one, (3, controlling the triplet correlations. 
Though the overlap of our best variational WF with the 
exact GS is rather poor, the average sign (s) is in general 
very much improved by the triplet term. Our interpre- 
tation is that short range many body correlations are 
very important to reproduce the relative phases of the 
GS on each Ising configuration. The optimal parameters 
for our initial guess ipy of the GS ipo are expected to 
be very weakly size-dependent but they are very difficult 
to determine accurately for large sizes. For a = 1 and 
N = 36, where ED is still possible, our best guess for the 
GS WF - with the maximum overlap and average sign 
- is slightly different from the one determined with the 
optimization of the energy. Since the forthcoming calcu- 
lations, which significantly improve the VMC, are more 
sensitive to the accuracy of the WF rather than to the 
one of the GS energy, henceforth we have chosen to work 
with j3 = 0.23 for all the system sizes. 
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FIG. 1. GS energy per site eo = Eo/N, in unit of J, 
as a function of the system size, obtained with VMC (full 
triangles), FN (empty dots) and SR with p = 7 (full dots) 
techniques. SW size scaling is assumed and short-dashed 
lines are linear fits against ljN 3 ^ 2 . The long-dashed line is 
the SW prediction, the empty triangle is the N — 36 ED 
result and the empty squares are data taken from Ref. Jl(]]. 

One way to get more accurate GS properties is to use 
the Green Function MC technique (GFMC). As in the 
fermionic case, for frustrated spin systems this numeri- 
cal method is plagued by the well-known sign problem. 
Recently, to alleviate the above mentioned instability, 
the Fixed-Node (FN) GFMC scheme @ has been in- 
troduced as a variational technique, typically much bet- 
ter than the conventional VMC. As shown in Fig. [|, and 



also pointed in Ref. for frustrated spin systems, this 
technique does not represent a significative advance com- 
pared to VMC, leading therefore to results biased by the 
variational ansatz. 
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FIG. 2. Short range spin correlation functions generated 
by H (a,b) and H 2 (c-g). 

In order to overcome this difficulty we have used a 
recently developed technique: GFMC with Stochastic 
Reconfiguration (SR) pH , which allows to release the 
FN approximation, in a controlled but approximate way, 
yielding, as shown in Fig. |l| much lower energies, even 
for the largest sizes where ED is not possible. In the ap- 
propriate limit pl| ] of large number of walkers and high 
frequency of SR, the residual bias introduced by the SR 
depends only on the number p of operators used to con- 
strain the GFMC Markov process. These constraints, 
analogously to the FN one, allow simulations without 
numerical instabilities. In principle the exact answer can 
be obtained, within statistical errors, provided p equals 
the huge Hilbert space dimension. Practically it is nec- 
essary to work with small p and an accurate selection of 
physically relevant operators is crucial. As can be eas- 
ily expected, the short range correlation functions S?Sj 

and (S^ S~ +S~ Sj~) contained in the Hamiltonian give 
a sizable improvement of the FN GS energy when they 
are put in the SR procedure. In order to be systematic 
we have included in the SR the short range correlations 
generated by H 2 (see Fig. |J), averaged over all spatial 
symmetries commuting with the Hamiltonian. This lo- 
cal correlations are particularly important to obtain quite 
accurate and reliable estimates not only of the GS en- 
ergy but also of the mixed average J2^] of the total spin 
square S? ot and of the order parameter w) 2 (defined as 
in Ref. |6j). These quantities are easily estimated within 
the GFMC technique and compared with the exact val- 
ues computed by ED for N = 36 in Tab. ||. In particular 
it is interesting that, starting from a variational WF with 
no definite spin, the GS singlet is systematically recov- 
ered by means of the SR technique. Furthermore, as it 
is shown in Fig. [l], the quality of our results is similar to 
the variational one obtained by P. Sindzingre et al. [ [To| , 
using a long range ordered RVB wavefunction. The latter 
approach is almost exact for small lattices, but the sign- 
problem is already present at the variational level, and 
the calculation has not been extended to high statistical 
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accuracy or to TV > 48. 

Having obtained an estimate for the GS energy, at least 
an order of magnitude more accurate than our best vari- 
ational guess, it appears possible to obtain physical fea- 
tures, such as a gap in the spin spectrum, that are not 
present at the variational level. For instance in the frus- 
trated J\ — Ji spin model, with the same technique and a 
similar accuracy, a gap in the spin spectrum was found in 
the thermodynamic limit, starting with a similar ordered 
and therefore gapless variational WF [|T) . 
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FIG. 3. Size scaling of the spin gap to the S = 3 excitation 
obtained with VMC, FN and SR (p = 7) techniques. The 
long-dashed line is the linear SW prediction and the solid line 
is the least-squares fit of the SR data for TV > 36. 

In the isotropic triangular AF, the gap to the first spin 
excitation is rather small. Furthermore, for the particu- 
lar choice of the guiding WF (||), the translational sym- 
metry of the Hamiltonian is preserved only if projected 
onto subspaces with total $f ot multiple of three. Such 
an 5 = 3 excitation belongs to the low-lying states of 
energy E$ and spin 5 of the ordered quantum AF, be- 
having as E s - E oc 5(5 + 1)/TV § . If instead E s - E 
remains finite for 5 = 3 and TV — > oo, this implies a dis- 
ordered GS. For all the above reasons we have studied 
the gap to the spin 5 = 3 excitation as a function of 
the system size. As it is shown in Fig. ||, for the lattice 
sizes for which a comparison with ED data is possible, 
the spin gap estimated with the SR technique is nearly 
exact. The importance to extend the numerical investi- 
gation to clusters large enough to allow a more reliable 
extrapolation is particularly evident in the same figure in 
which the TV = 12 and 36 exact data extrapolate linearly 
to a large finite value. This behavior, is certainly a finite 
size effect and it is corrected by the SR data for TV > 48, 
suggesting, strongly, a gapless excitation spectrum. 

As we have seen GFMC allows to obtain a very high 
statistical accuracy on the GS energy, but does not al- 
low to compute directly GS expectation values (4>o\0\ipo} 
p2| . A straightforward way is to perturb the Hamilto- 
nian with a term —AO , calculate the energy E(X) in 
presence of the perturbation and, by Hcllmann-Feynman 



theorem, estimate (ipo\0\^/ J o) = — dE(X)/dX\\=o with 
few computations at different small A's. A further 
complication for non exact calculations like the FN or 
SR, is that if the off-diagonal matrix elements O x '. x of 
the operator O (in the chosen basis) have the opposite 
sign of the product ijjy(x')ipv(x), they cannot be han- 
dled exactly within FN because these matrix elements 
change the nodes of ipv- A way to circumvent this dif- 
ficulty if to split the operator O in three contributions: 
O = D + + + 0~, where + (0~) is the operator with 
the same off-diagonal matrix elements of O when they 
have the same (opposite) signs of ipv{x')ipy(x), and zero 
otherwise, whereas D is the diagonal part of O. Then we 
can add to the Hamiltonian a contribution that does not 
change the nodes: H(X) = H — X(D + 20+) for A > 
and H(X) = H - X(D + 2 0~) for A < 0. It is easy to 
show that lim(E(-X) - E(X))/2X = (ip \6\ipo). 
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FIG. 4. Size scaling of the order parameter: VMC (full 
triangles), FN (empty dots), SR (full dots), exact data (empty 
triangles) and finite size linear SW (empty squares). The 
inset displays the A — » extrapolation for TV > 12. Lines are 
quadratic fits in all the plots. 

We plot in Fig. [Ito^ 2 estimated with this method using 
the FN and SR techniques. For the order parameter the 
inclusion of many short range correlations in the SR is 
not very important (see Tab. ||). Then, in order to min- 
imize the numerical effort, we have chosen to put in the 
SR conditions the first four correlation functions shown 
in Fig. ||, the order parameter itself and 5j ot . While the 
FN data extrapolate to a value not much lower than the 
variational result, the SR calculation provides a much 
more reliable estimate of the order parameter with no 
apparent loss of accuracy with increasing sizes. In this 
way we obtain for rrv a value well below the linear and 
the second order (which has actually a positive correction 
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fyl) SW predictions. This is partially in agreement with 
the conclusions of the finite temperature calculations |t]] 
suggesting a GS with a small but nonzero long range AF 
order and with series expansions Q indicating the tri- 
angular antiferromagnetic Heisenberg model to be likely 
ordered but close to a critical point. However in our simu- 
lation, which to our knowledge represents a first attempt 
to perform a systematic finite size scaling analysis of the 
order parameter, the value of rrv remains sizable and fi- 
nite, consistent with a gapless spectrum. This features 
could be also verified experimentally on the K/Si(lll):B 
interface (2^] which has turned out recently to be the first 
realization of a really bidimensional triangular AF. 

Though there is classical long range order, both the 
VMC and the SR approach show the crucial role of GS 
correlations defined on the smallest four spin clusters: in 
the variational calculation they are important to deter- 
mine the correct relative phases of the GS WF whereas 
in the latter more accurate approach this correlations al- 
low to obtain very accurate results for the energy and the 
spin gap and to restore the spin rotational invariance of 
the finite size GS. 
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TABLE I. Average sign, overlap, GS energy and its per- 
centage error obtained with the variational WF of Eq. ^ for 
N = 36 and some values of the easy-plane anisotropy a. The 
calculations were performed by summing exactly over all the 
configurations. 
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TABLE II. Variational estimate (VMC) and mixed aver- 
ages (FN, SR and Exact) of the GS energy per site, of the 
total spin square and of the AF order parameter. SR data 
are obtained using the first two (p = 2) , four (p — 4) and all 
(p = 7) the correlation functions shown in Fig. 
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